"""
Example 1.1

phi_xx + phi_yy = 0

phi = 0         on x = 0
phi = cos(pi*y) on x = 1
dphi/dn = 0 on y = 0 and y = 1 for 0 < x < 1.

Exact solution:

phi(x,y) = sinh(pi*x)*cos(pi*y)/sinh(pi)

"""

from numpy import sinh, cos, pi, zeros, exp
from bem import TwoDimensionalLaplaceBEMSolver

NO = 20
N = 4 * NO

dL = 1.0 / float(NO)

xb = zeros(N + 1)
yb = zeros(N + 1)

for i in range(NO):
    xb[i] = float(i) * dL
    yb[i] = 0.0

    xb[NO + i] = 1.0
    yb[NO + i] = xb[i]

    xb[2 * NO + i] = 1.0 - xb[i]
    yb[2 * NO + i] = 1.0

    xb[3 * NO + i] = 0.0
    yb[3 * NO + i] = 1.0 - xb[i]

    xb[N] = xb[0]
    yb[N] = yb[0]

bem = TwoDimensionalLaplaceBEMSolver(xb, yb)

for i in range(bem.number_of_boundary_elements):
    if i < NO:
        bem.bctypes[i] = 1
        bem.bcvalues[i] = 0.0
    elif i >= NO and i < 2 * NO:
        bem.bctypes[i] = 0
        bem.bcvalues[i] =  cos(pi * bem.ym[i])
    elif i >= 2 * NO and i < 3 * NO:
        bem.bctypes[i] = 1
        bem.bcvalues[i] = 0.0
    else:
        bem.bctypes[i] = 0
        bem.bcvalues[i] = 0.0

bem.solve()

def eval_bem_solution(x, y):
    return bem.evaluate_phi_at_interior_point(x, y)

def evaluate_exact_solution(x,y):
    return (sinh(pi * x) * cos(pi * y))/ sinh(pi)

if __name__ == "__main__":

    from numpy import zeros

    M = 50
    N = 50

    x = zeros((M,M),float)
    y = zeros((M,M),float)
    z = zeros((M,M),float)
    zz = zeros((M,M),float)
    xmin = 0.05
    xmax = 0.95
    dx = (xmax-xmin)/M
    ymin = 0.1
    ymax = 0.9
    dy = (ymax-ymin)/M

    for m in range(M):
        for n in range(N):
            x[m,n] = xmin + dx*n
            y[m,n] = ymin + dy*m
            z[m,n] = eval_bem_solution(x[m,n],y[m,n])

    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator, FormatStrFormatter
    import matplotlib.pyplot as plt


    plt.contourf(x, y, z)
    plt.show()
